I think the language could be cleaner, but this is how I have been processing zonal summaries. There are a couple steps to follow:

1. Open an FVCOM .nc file that has Lon/lat values to create a mesh. (this only needs to be done once, and will be buggy if its done each time since lat/lon sometimes are zeros on the THREDDS links)
2. Get geometric intersection of polygon and the mesh: get_poly_mesh_intersection()
3. Process the zonal average for all polygons in the intersection, this takes average of the three nodes: get_mesh_element_timeseries(), using get_average_from_nodes() on each element.
4. Get the zonal average by weighting each element-wise time series by its area (after intersection, so accounting for partial areas): get_area_weighted_zonal_average()

The function get_timeseries_for_poly() is a wrapper around doing steps 2-4.

Each of these need namespace bindings to be added, but I wanted to just contribute something to feel like I've done something.